library(readxl)
library(tidyverse)
library(plyr) # mapvalues()
library(descr) # freq()
library(DescTools) # statistic tests
library(modelsummary)
library(marginaleffects)
library(gridExtra) # combine graphs
library(lemon) # put legend in combined graphs
library(Hmisc) # correlation matrix & weighted sd
library(DescTools) # pseudo r-squared

PRRdataframe1 <- read_excel("Dataset.xlsx")

PRRdataframe2 <- na.omit(PRRdataframe1)

freq(PRRdataframe2$Party)

8546/2
4273/858 # 4.98
4273/7688 # 0.56

PRRdataframe2$df_weights <- mapvalues(PRRdataframe2$Party,
                                              from = c(0, 1),
                                              to = c(4.98, 0.56))

### PCAs on activism ###

# 1st run

PCA_dataframe <- PRRdataframe2[, c(16:20)]
pca_result <- prcomp(PCA_dataframe, scale = TRUE)
pca_result$rotation <- -1*pca_result$rotation
pca_result$x <- -1*pca_result$x

summary(pca_result)

# Scree plot
var_explained = pca_result$sdev^2 / sum(pca_result$sdev^2)

qplot(c(1:5), var_explained) + 
  geom_line() + 
  xlab("Principal Component") + 
  ylab("Variance Explained") +
  ggtitle("Scree Plot") +
  ylim(0, 1) +
  theme_bw()

# 2nd run

PCA_dataframe <- PRRdataframe2[, c(16, 18:20)]
pca_result <- prcomp(PCA_dataframe, scale = TRUE)
pca_result$rotation <- -1*pca_result$rotation
pca_result$x <- -1*pca_result$x

summary(pca_result)

# Scree plot
var_explained = pca_result$sdev^2 / sum(pca_result$sdev^2)

qplot(c(1:4), var_explained) + 
  geom_line() + 
  xlab("Principal Component") + 
  ylab("Variance Explained") +
  ggtitle("Scree plot") +
  ylim(0, 1) +
  theme_bw()

PCA_activism <- (pca_result$x[, 1]) # Extract scores of the first principal component (PC1)
PRRdataframe2$PCA_activism <- PCA_activism # Add the new variable to your original data frame


##############################
### DESCRIPTIVE STATISTICS ###
##############################

weightsss <- PRRdataframe1$party_wt

### AGGREGATE

freq(PRRdataframe1$meeting_num, plot = FALSE)
weighted.mean(PRRdataframe1$meeting_num, weightsss, na.rm=TRUE)
sqrt(wtd.var(PRRdataframe1$meeting_num, weightsss, na.rm=TRUE))
mean(PRRdataframe1$meeting_num, na.rm=TRUE)
sd(PRRdataframe1$meeting_num, na.rm=TRUE)

freq(PRRdataframe1$online_num, plot = FALSE)
weighted.mean(PRRdataframe1$online_num, weightsss, na.rm=TRUE)
sqrt(wtd.var(PRRdataframe1$online_num, weightsss, na.rm=TRUE))
mean(PRRdataframe1$online_num, na.rm=TRUE)
sd(PRRdataframe1$online_num, na.rm=TRUE)

freq(PRRdataframe1$flyers_num, plot = FALSE)
weighted.mean(PRRdataframe1$flyers_num, weightsss, na.rm=TRUE)
sqrt(wtd.var(PRRdataframe1$flyers_num, weightsss, na.rm=TRUE))
mean(PRRdataframe1$flyers_num, na.rm=TRUE)
sd(PRRdataframe1$flyers_num, na.rm=TRUE)

freq(PRRdataframe1$Partyoffice, plot = FALSE)
weighted.mean(PRRdataframe1$Partyoffice, weightsss, na.rm=TRUE)
sqrt(wtd.var(PRRdataframe1$Partyoffice, weightsss, na.rm=TRUE))
mean(PRRdataframe1$Partyoffice, na.rm=TRUE)
sd(PRRdataframe1$Partyoffice, na.rm=TRUE)

freq(PRRdataframe1$Candidature, plot = FALSE)
weighted.mean(PRRdataframe1$Candidature, weightsss, na.rm=TRUE)
sqrt(wtd.var(PRRdataframe1$Candidature, weightsss, na.rm=TRUE))
mean(PRRdataframe1$Candidature, na.rm=TRUE)
sd(PRRdataframe1$Candidature, na.rm=TRUE)

freq(PRRdataframe1$Gender, plot = FALSE)
weighted.mean(PRRdataframe1$Gender, weightsss, na.rm=TRUE)
sqrt(wtd.var(PRRdataframe1$Gender, weightsss, na.rm=TRUE))
mean(PRRdataframe1$Gender, na.rm=TRUE)
sd(PRRdataframe1$Gender, na.rm=TRUE)

freq(PRRdataframe1$age, plot = FALSE)
weighted.mean(PRRdataframe1$age, weightsss, na.rm=TRUE)
sqrt(wtd.var(PRRdataframe1$age, weightsss, na.rm=TRUE))
mean(PRRdataframe1$age, na.rm=TRUE)
sd(PRRdataframe1$age, na.rm=TRUE)

freq(PRRdataframe1$yearsmember, plot = FALSE)
weighted.mean(PRRdataframe1$yearsmember, weightsss, na.rm=TRUE)
sqrt(wtd.var(PRRdataframe1$yearsmember, weightsss, na.rm=TRUE))
mean(PRRdataframe1$yearsmember, na.rm=TRUE)
sd(PRRdataframe1$yearsmember, na.rm=TRUE)

freq(PRRdataframe1$Party, plot = FALSE)
weighted.mean(PRRdataframe1$Party, weightsss, na.rm=TRUE)
sqrt(wtd.var(PRRdataframe1$Party, weightsss, na.rm=TRUE))
mean(PRRdataframe1$Party, na.rm=TRUE)
sd(PRRdataframe1$Party, na.rm=TRUE)

freq(PRRdataframe1$Tertiaryedu, plot = FALSE)
weighted.mean(PRRdataframe1$Tertiaryedu, weightsss, na.rm=TRUE)
sqrt(wtd.var(PRRdataframe1$Tertiaryedu, weightsss, na.rm=TRUE))
mean(PRRdataframe1$Tertiaryedu, na.rm=TRUE)
sd(PRRdataframe1$Tertiaryedu, na.rm=TRUE)

freq(PRRdataframe1$Labourforce, plot = FALSE)
weighted.mean(PRRdataframe1$Labourforce, weightsss, na.rm=TRUE)
sqrt(wtd.var(PRRdataframe1$Labourforce, weightsss, na.rm=TRUE))
mean(PRRdataframe1$Labourforce, na.rm=TRUE)
sd(PRRdataframe1$Labourforce, na.rm=TRUE)

freq(PRRdataframe1$Policyinfluence, plot = FALSE)
weighted.mean(PRRdataframe1$Policyinfluence, weightsss, na.rm=TRUE)
sqrt(wtd.var(PRRdataframe1$Policyinfluence, weightsss, na.rm=TRUE))
mean(PRRdataframe1$Policyinfluence, na.rm=TRUE)
sd(PRRdataframe1$Policyinfluence, na.rm=TRUE)

freq(PRRdataframe1$Ambitioncand, plot = FALSE)
weighted.mean(PRRdataframe1$Ambitioncand, weightsss, na.rm=TRUE)
sqrt(wtd.var(PRRdataframe1$Ambitioncand, weightsss, na.rm=TRUE))
mean(PRRdataframe1$Ambitioncand, na.rm=TRUE)
sd(PRRdataframe1$Ambitioncand, na.rm=TRUE)

freq(PRRdataframe1$Newfriends, plot = FALSE)
weighted.mean(PRRdataframe1$Newfriends, weightsss, na.rm=TRUE)
sqrt(wtd.var(PRRdataframe1$Newfriends, weightsss, na.rm=TRUE))
mean(PRRdataframe1$Newfriends, na.rm=TRUE)
sd(PRRdataframe1$Newfriends, na.rm=TRUE)

freq(PRRdataframe1$Partyrecruitment, plot = FALSE)
weighted.mean(PRRdataframe1$Partyrecruitment, weightsss, na.rm=TRUE)
sqrt(wtd.var(PRRdataframe1$Partyrecruitment, weightsss, na.rm=TRUE))
mean(PRRdataframe1$Partyrecruitment, na.rm=TRUE)
sd(PRRdataframe1$Partyrecruitment, na.rm=TRUE)

freq(PRRdataframe1$Stigma, plot = FALSE)
weighted.mean(PRRdataframe1$Stigma, weightsss, na.rm=TRUE)
sqrt(wtd.var(PRRdataframe1$Stigma, weightsss, na.rm=TRUE))
mean(PRRdataframe1$Stigma, na.rm=TRUE)
sd(PRRdataframe1$Stigma, na.rm=TRUE)

### WOMEN ###

womendataframe <- PRRdataframe1 %>%
  filter(Gender == 1)

freq(womendataframe$online_num, plot = FALSE)
mean(womendataframe$online_num, na.rm=TRUE)
sd(womendataframe$online_num, na.rm=TRUE)

freq(womendataframe$meeting_num, plot = FALSE)
mean(womendataframe$meeting_num, na.rm=TRUE)
sd(womendataframe$meeting_num, na.rm=TRUE)

freq(womendataframe$flyers_num, plot = FALSE)
mean(womendataframe$flyers_num, na.rm=TRUE)
sd(womendataframe$flyers_num, na.rm=TRUE)

freq(womendataframe$Partyoffice, plot = FALSE)
mean(womendataframe$Partyoffice, na.rm=TRUE)
sd(womendataframe$Partyoffice, na.rm=TRUE)

freq(womendataframe$Candidature, plot = FALSE)
mean(womendataframe$Candidature, na.rm=TRUE)
sd(womendataframe$Candidature, na.rm=TRUE)

freq(womendataframe$age, plot = FALSE)
mean(womendataframe$age, na.rm=TRUE)
sd(womendataframe$age, na.rm=TRUE)

freq(womendataframe$Tertiaryedu, plot = FALSE)
mean(womendataframe$Tertiaryedu, na.rm=TRUE)
sd(womendataframe$Tertiaryedu, na.rm=TRUE)

freq(womendataframe$Labourforce, plot = FALSE)
mean(womendataframe$Labourforce, na.rm=TRUE)
sd(womendataframe$Labourforce, na.rm=TRUE)

freq(womendataframe$yearsmember, plot = FALSE)
mean(womendataframe$yearsmember, na.rm=TRUE)
sd(womendataframe$yearsmember, na.rm=TRUE)

freq(womendataframe$Party, plot = FALSE)
mean(womendataframe$Party, na.rm=TRUE)
sd(womendataframe$Party, na.rm=TRUE)

freq(womendataframe$Policyinfluence, plot = FALSE)
mean(womendataframe$Policyinfluence, na.rm=TRUE)
sd(womendataframe$Policyinfluence, na.rm=TRUE)

freq(womendataframe$Ambitioncand, plot = FALSE)
mean(womendataframe$Ambitioncand, na.rm=TRUE)
sd(womendataframe$Ambitioncand, na.rm=TRUE)

freq(womendataframe$Partyrecruitment, plot = FALSE)
mean(womendataframe$Partyrecruitment, na.rm=TRUE)
sd(womendataframe$Partyrecruitment, na.rm=TRUE)

freq(womendataframe$Newfriends, plot = FALSE)
mean(womendataframe$Newfriends, na.rm=TRUE)
sd(womendataframe$Newfriends, na.rm=TRUE)

freq(womendataframe$Stigma, plot = FALSE)
mean(womendataframe$Stigma, na.rm=TRUE)
sd(womendataframe$Stigma, na.rm=TRUE)

### MEN ###

mendataframe <- PRRdataframe1 %>%
  filter(Gender == 0)

freq(mendataframe$online_num, plot = FALSE)
mean(mendataframe$online_num, na.rm=TRUE)
sd(mendataframe$online_num, na.rm=TRUE)

freq(mendataframe$meeting_num, plot = FALSE)
mean(mendataframe$meeting_num, na.rm=TRUE)
sd(mendataframe$meeting_num, na.rm=TRUE)

freq(mendataframe$flyers_num, plot = FALSE)
mean(mendataframe$flyers_num, na.rm=TRUE)
sd(mendataframe$flyers_num, na.rm=TRUE)

freq(mendataframe$Partyoffice, plot = FALSE)
mean(mendataframe$Partyoffice, na.rm=TRUE)
sd(mendataframe$Partyoffice, na.rm=TRUE)

freq(mendataframe$age, plot = FALSE)
mean(mendataframe$age, na.rm=TRUE)
sd(mendataframe$age, na.rm=TRUE)

freq(mendataframe$Tertiaryedu, plot = FALSE)
mean(mendataframe$Tertiaryedu, na.rm=TRUE)
sd(mendataframe$Tertiaryedu, na.rm=TRUE)

freq(mendataframe$Labourforce, plot = FALSE)
mean(mendataframe$Labourforce, na.rm=TRUE)
sd(mendataframe$Labourforce, na.rm=TRUE)

freq(mendataframe$yearsmember, plot = FALSE)
mean(mendataframe$yearsmember, na.rm=TRUE)
sd(mendataframe$yearsmember, na.rm=TRUE)

freq(mendataframe$Party, plot = FALSE)
mean(mendataframe$Party, na.rm=TRUE)
sd(mendataframe$Party, na.rm=TRUE)

freq(mendataframe$Policyinfluence, plot = FALSE)
mean(mendataframe$Policyinfluence, na.rm=TRUE)
sd(mendataframe$Policyinfluence, na.rm=TRUE)

freq(mendataframe$Ambitioncand, plot = FALSE)
mean(mendataframe$Ambitioncand, na.rm=TRUE)
sd(mendataframe$Ambitioncand, na.rm=TRUE)

freq(mendataframe$Partyrecruitment, plot = FALSE)
mean(mendataframe$Partyrecruitment, na.rm=TRUE)
sd(mendataframe$Partyrecruitment, na.rm=TRUE)

freq(mendataframe$Newfriends, plot = FALSE)
mean(mendataframe$Newfriends, na.rm=TRUE)
sd(mendataframe$Newfriends, na.rm=TRUE)

freq(mendataframe$Stigma, plot = FALSE)
mean(mendataframe$Stigma, na.rm=TRUE)
sd(mendataframe$Stigma, na.rm=TRUE)

### LEAGUE ###

League_data <- PRRdataframe1 %>%
  filter(Party == 0)

freq(League_data$online_num, plot = FALSE)
mean(League_data$online_num, na.rm=TRUE)
sd(League_data$online_num, na.rm=TRUE)

freq(League_data$meeting_num, plot = FALSE)
mean(League_data$meeting_num, na.rm=TRUE)
sd(League_data$meeting_num, na.rm=TRUE)

freq(League_data$flyers_num, plot = FALSE)
mean(League_data$flyers_num, na.rm=TRUE)
sd(League_data$flyers_num, na.rm=TRUE)

freq(League_data$Partyoffice, plot = FALSE)
mean(League_data$Partyoffice, na.rm=TRUE)
sd(League_data$Partyoffice, na.rm=TRUE)

freq(League_data$Candidature, plot = FALSE)
mean(League_data$Candidature, na.rm=TRUE)
sd(League_data$Candidature, na.rm=TRUE)

freq(League_data$Gender, plot = FALSE)
mean(League_data$Gender, na.rm=TRUE)
sd(League_data$Gender, na.rm=TRUE)

freq(League_data$age, plot = FALSE)
mean(League_data$age, na.rm=TRUE)
sd(League_data$age, na.rm=TRUE)

freq(League_data$Tertiaryedu, plot = FALSE)
mean(League_data$Tertiaryedu, na.rm=TRUE)
sd(League_data$Tertiaryedu, na.rm=TRUE)

freq(League_data$Labourforce, plot = FALSE)
mean(League_data$Labourforce, na.rm=TRUE)
sd(League_data$Labourforce, na.rm=TRUE)

freq(League_data$yearsmember, plot = FALSE)
mean(League_data$yearsmember, na.rm=TRUE)
sd(League_data$yearsmember, na.rm=TRUE)

freq(League_data$Policyinfluence, plot = FALSE)
mean(League_data$Policyinfluence, na.rm=TRUE)
sd(League_data$Policyinfluence, na.rm=TRUE)

freq(League_data$Ambitioncand, plot = FALSE)
mean(League_data$Ambitioncand, na.rm=TRUE)
sd(League_data$Ambitioncand, na.rm=TRUE)

freq(League_data$Partyrecruitment, plot = FALSE)
mean(League_data$Partyrecruitment, na.rm=TRUE)
sd(League_data$Partyrecruitment, na.rm=TRUE)

freq(League_data$Newfriends, plot = FALSE)
mean(League_data$Newfriends, na.rm=TRUE)
sd(League_data$Newfriends, na.rm=TRUE)

freq(League_data$Stigma, plot = FALSE)
mean(League_data$Stigma, na.rm=TRUE)
sd(League_data$Stigma, na.rm=TRUE)

### SD ###

SD_data <- PRRdataframe1 %>%
  filter(Party == 1)

freq(SD_data$online_num, plot = FALSE)
mean(SD_data$online_num, na.rm=TRUE)
sd(SD_data$online_num, na.rm=TRUE)

freq(SD_data$meeting_num, plot = FALSE)
mean(SD_data$meeting_num, na.rm=TRUE)
sd(SD_data$meeting_num, na.rm=TRUE)

freq(SD_data$flyers_num, plot = FALSE)
mean(SD_data$flyers_num, na.rm=TRUE)
sd(SD_data$flyers_num, na.rm=TRUE)

freq(SD_data$Partyoffice, plot = FALSE)
mean(SD_data$Partyoffice, na.rm=TRUE)
sd(SD_data$Partyoffice, na.rm=TRUE)

freq(SD_data$Candidature, plot = FALSE)
mean(SD_data$Candidature, na.rm=TRUE)
sd(SD_data$Candidature, na.rm=TRUE)

freq(SD_data$Gender, plot = FALSE)
mean(SD_data$Gender, na.rm=TRUE)
sd(SD_data$Gender, na.rm=TRUE)

freq(SD_data$age, plot = FALSE)
mean(SD_data$age, na.rm=TRUE)
sd(SD_data$age, na.rm=TRUE)

freq(SD_data$Tertiaryedu, plot = FALSE)
mean(SD_data$Tertiaryedu, na.rm=TRUE)
sd(SD_data$Tertiaryedu, na.rm=TRUE)

freq(SD_data$Labourforce, plot = FALSE)
mean(SD_data$Labourforce, na.rm=TRUE)
sd(SD_data$Labourforce, na.rm=TRUE)

freq(SD_data$yearsmember, plot = FALSE)
mean(SD_data$yearsmember, na.rm=TRUE)
sd(SD_data$yearsmember, na.rm=TRUE)

freq(SD_data$Policyinfluence, plot = FALSE)
mean(SD_data$Policyinfluence, na.rm=TRUE)
sd(SD_data$Policyinfluence, na.rm=TRUE)

freq(SD_data$Ambitioncand, plot = FALSE)
mean(SD_data$Ambitioncand, na.rm=TRUE)
sd(SD_data$Ambitioncand, na.rm=TRUE)

freq(SD_data$Partyrecruitment, plot = FALSE)
mean(SD_data$Partyrecruitment, na.rm=TRUE)
sd(SD_data$Partyrecruitment, na.rm=TRUE)

freq(SD_data$Newfriends, plot = FALSE)
mean(SD_data$Newfriends, na.rm=TRUE)
sd(SD_data$Newfriends, na.rm=TRUE)

freq(SD_data$Stigma, plot = FALSE)
mean(SD_data$Stigma, na.rm=TRUE)
sd(SD_data$Stigma, na.rm=TRUE)

#########################
### MULTICOLLINEARITY ###
#########################

newdataframe <- PRRdataframe2[, c(2, 4, 6:11, 13:15)]

corrmatrix <- rcorr(as.matrix(newdataframe))
corrmatrix$r

corrmatrix <- as.data.frame(corrmatrix$r)


##############
### MODELS ###
##############

library(MASS)

tidy_custom.polr <- function(x, ...) {
  s <- lmtest::coeftest(x)
  out <- data.frame(
    term = row.names(s),
    p.value = s[, "Pr(>|t|)"])
  out
}

coefficientsnames <- c("Party" = "Party (SD)",
                       "yearsmember" = "Years of membership",
                       "Labourforce" = "Active in the workforce",
                       "Tertiaryedu" = "Tertiary education",
                       "agesquared" = "Age squared",
                       "age" = "Age",
                       "Stigma" = "Stigma",
                       "Newfriends" = "New friends since joining",
                       "Partyrecruitment" = "Recruited by the party",
                       "Ambitioncand" = "Electoral ambition",
                       "Policyinfluence" = "Political efficacy",
                       "Gender" = "Gender (Woman)")


### STEPWISE MODELS ###

# Meetings

model_1 <- polr(as.ordered(meeting_num) ~ Gender,
                data = PRRdataframe2, weights = df_weights)
modelsummary(model_1, stars = TRUE, vcov = ~Party, coef_rename = coefficientsnames)

model_2 <- polr(as.ordered(meeting_num) ~ Gender + Policyinfluence + Ambitioncand,
                data = PRRdataframe2, weights = df_weights)
modelsummary(model_2, stars = TRUE, vcov = ~Party, coef_rename = coefficientsnames)

model_3 <- polr(as.ordered(meeting_num) ~ Gender + Policyinfluence + Ambitioncand +
                  Partyrecruitment + Newfriends,
                data = PRRdataframe2, weights = df_weights)
modelsummary(model_3, stars = TRUE, vcov = ~Party, coef_rename = coefficientsnames)

model_4 <- polr(as.ordered(meeting_num) ~ Gender + Policyinfluence + Ambitioncand +
                  Partyrecruitment + Newfriends + Stigma,
                data = PRRdataframe2, weights = df_weights)
modelsummary(model_4, stars = TRUE, vcov = ~Party, coef_rename = coefficientsnames)

model_5 <- polr(as.ordered(meeting_num) ~ Gender + Policyinfluence + Ambitioncand +
                  Partyrecruitment + Newfriends + Stigma +
                  age + agesquared + Tertiaryedu + Labourforce + yearsmember + Party,
                data = PRRdataframe2, weights = df_weights)
modelsummary(model_5, stars = TRUE, vcov = ~Party, coef_rename = coefficientsnames)

meeting_table <- list(model_1, model_2, model_3, model_4, model_5)

modelsummary(meeting_table, stars = TRUE, vcov = ~Party, coef_rename = coefficientsnames,
             output = "Attending meetings.docx")

# Online

model_7 <- polr(as.ordered(online_num) ~ Gender,
                data = PRRdataframe2, weights = df_weights)
modelsummary(model_7, stars = TRUE, vcov = ~Party, coef_rename = coefficientsnames)

model_8 <- polr(as.ordered(online_num) ~ Gender + Policyinfluence + Ambitioncand,
                data = PRRdataframe2, weights = df_weights)
modelsummary(model_8, stars = TRUE, vcov = ~Party, coef_rename = coefficientsnames)

model_9 <- polr(as.ordered(online_num) ~ Gender + Policyinfluence + Ambitioncand +
                  Partyrecruitment + Newfriends,
                data = PRRdataframe2, weights = df_weights)
modelsummary(model_9, stars = TRUE, vcov = ~Party, coef_rename = coefficientsnames)

model_10 <- polr(as.ordered(online_num) ~ Gender + Policyinfluence + Ambitioncand +
                   Partyrecruitment + Newfriends + Stigma,
                 data = PRRdataframe2, weights = df_weights)
modelsummary(model_10, stars = TRUE, vcov = ~Party, coef_rename = coefficientsnames)

model_11 <- polr(as.ordered(online_num) ~ Gender + Policyinfluence + Ambitioncand +
                   Partyrecruitment + Newfriends + Stigma +
                   age + agesquared + Tertiaryedu + Labourforce + yearsmember + Party,
                 data = PRRdataframe2, weights = df_weights)
modelsummary(model_11, stars = TRUE, vcov = ~Party, coef_rename = coefficientsnames)

online_table <- list(model_7, model_8, model_9, model_10, model_11)

modelsummary(online_table, stars = TRUE, vcov = ~Party, coef_rename = coefficientsnames,
             output = "Online.docx")

# Flyers

model_13 <- polr(as.ordered(flyers_num) ~ Gender,
                 data = PRRdataframe2, weights = df_weights)
modelsummary(model_13, stars = TRUE, vcov = ~Party, coef_rename = coefficientsnames)

model_14 <- polr(as.ordered(flyers_num) ~ Gender + Policyinfluence + Ambitioncand,
                 data = PRRdataframe2, weights = df_weights)
modelsummary(model_14, stars = TRUE, vcov = ~Party, coef_rename = coefficientsnames)

model_15 <- polr(as.ordered(flyers_num) ~ Gender + Policyinfluence + Ambitioncand +
                   Partyrecruitment + Newfriends,
                 data = PRRdataframe2, weights = df_weights)
modelsummary(model_15, stars = TRUE, vcov = ~Party, coef_rename = coefficientsnames)

model_16 <- polr(as.ordered(flyers_num) ~ Gender + Policyinfluence + Ambitioncand +
                   Partyrecruitment + Newfriends + Stigma,
                 data = PRRdataframe2, weights = df_weights)
modelsummary(model_16, stars = TRUE, vcov = ~Party, coef_rename = coefficientsnames)

model_17 <- polr(as.ordered(flyers_num) ~ Gender + Policyinfluence + Ambitioncand +
                   Partyrecruitment + Newfriends + Stigma +
                   age + agesquared + Tertiaryedu + Labourforce + yearsmember + Party,
                 data = PRRdataframe2, weights = df_weights)
modelsummary(model_17, stars = TRUE, vcov = ~Party, coef_rename = coefficientsnames)

flyers_table <- list(model_13, model_14, model_15, model_16, model_17)

modelsummary(flyers_table, stars = TRUE, vcov = ~Party, coef_rename = coefficientsnames,
             output = "Flyers.docx")

# Party office

model_19 <- glm(Partyoffice ~ Gender,
                data = PRRdataframe2, weights = df_weights, family = "binomial")
modelsummary(model_19, stars = TRUE, vcov = ~Party, coef_rename = coefficientsnames)
PseudoR2(model_19)

model_20 <- glm(Partyoffice ~ Gender + Policyinfluence + Ambitioncand,
                data = PRRdataframe2, weights = df_weights, family = "binomial")
modelsummary(model_20, stars = TRUE, vcov = ~Party, coef_rename = coefficientsnames)
PseudoR2(model_20)

model_21 <- glm(Partyoffice ~ Gender + Policyinfluence + Ambitioncand +
                  Partyrecruitment + Newfriends,
                data = PRRdataframe2, weights = df_weights, family = "binomial")
modelsummary(model_21, stars = TRUE, vcov = ~Party, coef_rename = coefficientsnames)
PseudoR2(model_21)

model_22 <- glm(Partyoffice ~ Gender + Policyinfluence + Ambitioncand +
                  Partyrecruitment + Newfriends + Stigma,
                data = PRRdataframe2, weights = df_weights, family = "binomial")
modelsummary(model_22, stars = TRUE, vcov = ~Party, coef_rename = coefficientsnames)
PseudoR2(model_22)

model_23 <- glm(Partyoffice ~ Gender + Policyinfluence + Ambitioncand +
                  Partyrecruitment + Newfriends + Stigma +
                  age + agesquared + Tertiaryedu + Labourforce + yearsmember + Party,
                data = PRRdataframe2, weights = df_weights, family = "binomial")
modelsummary(model_23, stars = TRUE, vcov = ~Party, coef_rename = coefficientsnames)
PseudoR2(model_23)

partyoffice_table <- list(model_19, model_20, model_21, model_22, model_23)

modelsummary(partyoffice_table, stars = TRUE, vcov = ~Party, coef_rename = coefficientsnames,
             output = "Party office.docx")

# Candidature

model_25 <- glm(Candidature ~ Gender,
                data = PRRdataframe2, weights = df_weights, family = "binomial")
modelsummary(model_25, stars = TRUE, vcov = ~Party, coef_rename = coefficientsnames)
PseudoR2(model_25)

model_26 <- glm(Candidature ~ Gender + Policyinfluence + Ambitioncand,
                data = PRRdataframe2, weights = df_weights, family = "binomial")
modelsummary(model_26, stars = TRUE, vcov = ~Party, coef_rename = coefficientsnames)
PseudoR2(model_26)

model_27 <- glm(Candidature ~ Gender + Policyinfluence + Ambitioncand +
                  Partyrecruitment + Newfriends,
                data = PRRdataframe2, weights = df_weights, family = "binomial")
modelsummary(model_27, stars = TRUE, vcov = ~Party, coef_rename = coefficientsnames)
PseudoR2(model_27)

model_28 <- glm(Candidature ~ Gender + Policyinfluence + Ambitioncand +
                  Partyrecruitment + Newfriends + Stigma,
                data = PRRdataframe2, weights = df_weights, family = "binomial")
modelsummary(model_28, stars = TRUE, vcov = ~Party, coef_rename = coefficientsnames)
PseudoR2(model_28)

model_29 <- glm(Candidature ~ Gender + Policyinfluence + Ambitioncand +
                  Partyrecruitment + Newfriends + Stigma +
                  age + agesquared + Tertiaryedu + Labourforce + yearsmember + Party,
                data = PRRdataframe2, weights = df_weights, family = "binomial")
modelsummary(model_29, stars = TRUE, vcov = ~Party, coef_rename = coefficientsnames)
PseudoR2(model_29)

candidature_table <- list(model_25, model_26, model_27, model_28, model_29)

modelsummary(candidature_table, stars = TRUE, vcov = ~Party, coef_rename = coefficientsnames,
             output = "Candidature.docx")

### BIVARIATE ###

bivariates_table <- list(model_7, model_1, model_13, model_19, model_25)

modelsummary(bivariates_table, stars = TRUE, vcov = ~Party, coef_rename = coefficientsnames,
             output = "Bivariate regressions.docx")


#########################
### ROBUSTNESS CHECKS ###
#########################

### ONLY LEAGUE ###

Leaguedataframe <- PRRdataframe2 %>%
  filter(Party == 0)

# Full models

rc1 <- polr(as.ordered(meeting_num) ~ Gender + Policyinfluence + Ambitioncand +
              Partyrecruitment + Newfriends + Stigma +
              age + agesquared + Tertiaryedu + Labourforce + yearsmember,
               data = Leaguedataframe)
modelsummary(rc1, stars = TRUE, vcov = ~county, coef_rename = coefficientsnames)

rc2 <- polr(as.ordered(online_num) ~ Gender + Policyinfluence + Ambitioncand +
              Partyrecruitment + Newfriends + Stigma +
              age + agesquared + Tertiaryedu + Labourforce + yearsmember,
               data = Leaguedataframe)
modelsummary(rc2, stars = TRUE, vcov = ~county, coef_rename = coefficientsnames)

rc3 <- polr(as.ordered(flyers_num) ~ Gender + Policyinfluence + Ambitioncand +
              Partyrecruitment + Newfriends + Stigma +
              age + agesquared + Tertiaryedu + Labourforce + yearsmember,
               data = Leaguedataframe)
modelsummary(rc3, stars = TRUE, vcov = ~county, coef_rename = coefficientsnames)

rc4 <- glm(Partyoffice ~ Gender + Policyinfluence + Ambitioncand +
             Partyrecruitment + Newfriends + Stigma +
             age + agesquared + Tertiaryedu + Labourforce + yearsmember,
              data = Leaguedataframe, family = "binomial")
modelsummary(rc4, stars = TRUE, vcov = ~county, coef_rename = coefficientsnames)
PseudoR2(rc4)

rc5 <- glm(Candidature ~ Gender + Policyinfluence + Ambitioncand +
             Partyrecruitment + Newfriends + Stigma +
             age + agesquared + Tertiaryedu + Labourforce + yearsmember,
              data = Leaguedataframe, family = "binomial")
modelsummary(rc5, stars = TRUE, vcov = ~county, coef_rename = coefficientsnames)
PseudoR2(rc5)

rctable1 <- list(rc2, rc1, rc3, rc4, rc5)

modelsummary(rctable1, stars = TRUE, vcov = ~county, coef_rename = coefficientsnames,
             output = "Table RC Lega.docx")


### ONLY SD ###

SDdataframe <- PRRdataframe2 %>%
  filter(Party == 1)

rc6 <- polr(as.ordered(meeting_num) ~ Gender + Policyinfluence + Ambitioncand +
              Partyrecruitment + Newfriends + Stigma +
              age + agesquared + Tertiaryedu + Labourforce + yearsmember,
            data = SDdataframe)
modelsummary(rc6, stars = TRUE, vcov = ~county, coef_rename = coefficientsnames)

rc7 <- polr(as.ordered(online_num) ~ Gender + Policyinfluence + Ambitioncand +
              Partyrecruitment + Newfriends + Stigma +
              age + agesquared + Tertiaryedu + Labourforce + yearsmember,
            data = SDdataframe)
modelsummary(rc7, stars = TRUE, vcov = ~county, coef_rename = coefficientsnames)

rc8 <- polr(as.ordered(flyers_num) ~ Gender + Policyinfluence + Ambitioncand +
              Partyrecruitment + Newfriends + Stigma +
              age + agesquared + Tertiaryedu + Labourforce + yearsmember,
            data = SDdataframe)
modelsummary(rc8, stars = TRUE, vcov = ~county, coef_rename = coefficientsnames)

rc9 <- glm(Partyoffice ~ Gender + Policyinfluence + Ambitioncand +
             Partyrecruitment + Newfriends + Stigma +
             age + agesquared + Tertiaryedu + Labourforce + yearsmember,
           data = SDdataframe, family = "binomial")
modelsummary(rc9, stars = TRUE, vcov = ~county, coef_rename = coefficientsnames)
PseudoR2(rc9)

rc10 <- glm(Candidature ~ Gender + Policyinfluence + Ambitioncand +
              Partyrecruitment + Newfriends + Stigma +
              age + agesquared + Tertiaryedu + Labourforce + yearsmember,
           data = SDdataframe, family = "binomial")
modelsummary(rc10, stars = TRUE, vcov = ~county, coef_rename = coefficientsnames)
PseudoR2(rc10)

rctable2 <- list(rc7, rc6, rc8, rc9, rc10)

modelsummary(rctable2, stars = TRUE, vcov = ~county, coef_rename = coefficientsnames,
             output = "Table RC SD.docx")


### LEAGUE: ONLY MAILING LIST ###

Leaguedataframe$mailinglist <- mapvalues(Leaguedataframe$county,
                                from = c("Abruzzo", "Campania", "Emilia-Romagna", "Friuli-Venezia Giulia",
                                         "Lazio", "Liguria", "Lombardia", "Marche", "Molise", "Piemonte",
                                         "Puglia", "Sardegna", "Sicilia", "Toscana", "Umbria", "Veneto"),
                                to = c("No", "No", "Yes", "No", "No", "No", "No", "No", "No", "No", "No",
                                       "Yes", "Yes", "Yes", "Yes", "Yes"))

League_mail <- Leaguedataframe %>%
  filter(mailinglist == "Yes")

rc11 <- polr(as.ordered(meeting_num) ~ Gender + Policyinfluence + Ambitioncand +
               Partyrecruitment + Newfriends + Stigma +
               age + agesquared + Tertiaryedu + Labourforce + yearsmember,
            data = League_mail)
modelsummary(rc11, stars = TRUE, vcov = ~county, coef_rename = coefficientsnames)

rc12 <- polr(as.ordered(online_num) ~ Gender + Policyinfluence + Ambitioncand +
               Partyrecruitment + Newfriends + Stigma +
               age + agesquared + Tertiaryedu + Labourforce + yearsmember,
            data = League_mail)
modelsummary(rc12, stars = TRUE, vcov = ~county, coef_rename = coefficientsnames)

rc13 <- polr(as.ordered(flyers_num) ~ Gender + Policyinfluence + Ambitioncand +
               Partyrecruitment + Newfriends + Stigma +
               age + agesquared + Tertiaryedu + Labourforce + yearsmember,
            data = League_mail)
modelsummary(rc13, stars = TRUE, vcov = ~county, coef_rename = coefficientsnames)

rc14 <- glm(Partyoffice ~ Gender + Policyinfluence + Ambitioncand +
              Partyrecruitment + Newfriends + Stigma +
              age + agesquared + Tertiaryedu + Labourforce + yearsmember,
           data = League_mail, family = "binomial")
modelsummary(rc14, stars = TRUE, vcov = ~county, coef_rename = coefficientsnames)
PseudoR2(rc14)

rc15 <- glm(Candidature ~ Gender + Policyinfluence + Ambitioncand +
              Partyrecruitment + Newfriends + Stigma +
              age + agesquared + Tertiaryedu + Labourforce + yearsmember,
           data = League_mail, family = "binomial")
modelsummary(rc15, stars = TRUE, vcov = ~county, coef_rename = coefficientsnames)
PseudoR2(rc15)

rctable3 <- list(rc12, rc11, rc13, rc14, rc15)

modelsummary(rctable3, stars = TRUE, vcov = ~county, coef_rename = coefficientsnames,
             output = "Table RC Lega - Mailing list.docx")


### SD: NO SKANE ###

SDdfnoska <- SDdataframe %>%
  filter(county != "Skäne")

rc16 <- polr(as.ordered(meeting_num) ~ Gender + Policyinfluence + Ambitioncand +
               Partyrecruitment + Newfriends + Stigma +
               age + agesquared + Tertiaryedu + Labourforce + yearsmember,
             data = SDdfnoska)
modelsummary(rc16, stars = TRUE, vcov = ~county, coef_rename = coefficientsnames)

rc17 <- polr(as.ordered(online_num) ~ Gender + Policyinfluence + Ambitioncand +
               Partyrecruitment + Newfriends + Stigma +
               age + agesquared + Tertiaryedu + Labourforce + yearsmember,
             data = SDdfnoska)
modelsummary(rc17, stars = TRUE, vcov = ~county, coef_rename = coefficientsnames)

rc18 <- polr(as.ordered(flyers_num) ~ Gender + Policyinfluence + Ambitioncand +
               Partyrecruitment + Newfriends + Stigma +
               age + agesquared + Tertiaryedu + Labourforce + yearsmember,
             data = SDdfnoska)
modelsummary(rc18, stars = TRUE, vcov = ~county, coef_rename = coefficientsnames)

rc19 <- glm(Partyoffice ~ Gender + Policyinfluence + Ambitioncand +
              Partyrecruitment + Newfriends + Stigma +
              age + agesquared + Tertiaryedu + Labourforce + yearsmember,
            data = SDdfnoska, family = "binomial")
modelsummary(rc19, stars = TRUE, vcov = ~county, coef_rename = coefficientsnames)
PseudoR2(rc19)

rc20 <- glm(Candidature ~ Gender + Policyinfluence + Ambitioncand +
              Partyrecruitment + Newfriends + Stigma +
              age + agesquared + Tertiaryedu + Labourforce + yearsmember,
            data = SDdfnoska, family = "binomial")
modelsummary(rc20, stars = TRUE, vcov = ~county, coef_rename = coefficientsnames)
PseudoR2(rc20)

rctable4 <- list(rc17, rc16, rc18, rc19, rc20)

modelsummary(rctable4, stars = TRUE, vcov = ~county, coef_rename = coefficientsnames,
             output = "Table RC SD - No Skane.docx")


### SD: NO STOCKHOLM ###

SDdfnosto <- SDdataframe %>%
  filter(county != "Stockholm")

rc21 <- polr(as.ordered(meeting_num) ~ Gender + Policyinfluence + Ambitioncand +
               Partyrecruitment + Newfriends + Stigma +
               age + agesquared + Tertiaryedu + Labourforce + yearsmember,
             data = SDdfnosto)
modelsummary(rc21, stars = TRUE, vcov = ~county, coef_rename = coefficientsnames)

rc22 <- polr(as.ordered(online_num) ~ Gender + Policyinfluence + Ambitioncand +
               Partyrecruitment + Newfriends + Stigma +
               age + agesquared + Tertiaryedu + Labourforce + yearsmember,
             data = SDdfnosto)
modelsummary(rc22, stars = TRUE, vcov = ~county, coef_rename = coefficientsnames)

rc23 <- polr(as.ordered(flyers_num) ~ Gender + Policyinfluence + Ambitioncand +
               Partyrecruitment + Newfriends + Stigma +
               age + agesquared + Tertiaryedu + Labourforce + yearsmember,
             data = SDdfnosto)
modelsummary(rc23, stars = TRUE, vcov = ~county, coef_rename = coefficientsnames)

rc24 <- glm(Partyoffice ~ Gender + Policyinfluence + Ambitioncand +
              Partyrecruitment + Newfriends + Stigma +
              age + agesquared + Tertiaryedu + Labourforce + yearsmember,
            data = SDdfnosto, family = "binomial")
modelsummary(rc24, stars = TRUE, vcov = ~county, coef_rename = coefficientsnames)
PseudoR2(rc24)

rc25 <- glm(Candidature ~ Gender + Policyinfluence + Ambitioncand +
              Partyrecruitment + Newfriends + Stigma +
              age + agesquared + Tertiaryedu + Labourforce + yearsmember,
            data = SDdfnosto, family = "binomial")
modelsummary(rc25, stars = TRUE, vcov = ~county, coef_rename = coefficientsnames)
PseudoR2(rc25)

rctable5 <- list(rc22, rc21, rc23, rc24, rc25)

modelsummary(rctable5, stars = TRUE, vcov = ~county, coef_rename = coefficientsnames,
             output = "Table RC SD - No Stockholm.docx")

### SD: NO VG ###

SDdfnoVG <- SDdataframe %>%
  filter(county != "Västra Götaland")

rc26 <- polr(as.ordered(meeting_num) ~ Gender + Policyinfluence + Ambitioncand +
               Partyrecruitment + Newfriends + Stigma +
               age + agesquared + Tertiaryedu + Labourforce + yearsmember,
             data = SDdfnoVG)
modelsummary(rc26, stars = TRUE, vcov = ~county, coef_rename = coefficientsnames)

rc27 <- polr(as.ordered(online_num) ~ Gender + Policyinfluence + Ambitioncand +
               Partyrecruitment + Newfriends + Stigma +
               age + agesquared + Tertiaryedu + Labourforce + yearsmember,
             data = SDdfnoVG)
modelsummary(rc27, stars = TRUE, vcov = ~county, coef_rename = coefficientsnames)

rc28 <- polr(as.ordered(flyers_num) ~ Gender + Policyinfluence + Ambitioncand +
               Partyrecruitment + Newfriends + Stigma +
               age + agesquared + Tertiaryedu + Labourforce + yearsmember,
             data = SDdfnoVG)
modelsummary(rc28, stars = TRUE, vcov = ~county, coef_rename = coefficientsnames)

rc29 <- glm(Partyoffice ~ Gender + Policyinfluence + Ambitioncand +
              Partyrecruitment + Newfriends + Stigma +
              age + agesquared + Tertiaryedu + Labourforce + yearsmember,
            data = SDdfnoVG, family = "binomial")
modelsummary(rc29, stars = TRUE, vcov = ~county, coef_rename = coefficientsnames)
PseudoR2(rc29)

rc30 <- glm(Candidature ~ Gender + Policyinfluence + Ambitioncand +
              Partyrecruitment + Newfriends + Stigma +
              age + agesquared + Tertiaryedu + Labourforce + yearsmember,
            data = SDdfnoVG, family = "binomial")
modelsummary(rc30, stars = TRUE, vcov = ~county, coef_rename = coefficientsnames)
PseudoR2(rc30)

rctable6 <- list(rc27, rc26, rc28, rc29, rc30)

modelsummary(rctable6, stars = TRUE, vcov = ~county, coef_rename = coefficientsnames,
             output = "Table RC SD - No VG.docx")


### SD: NO 3 ###

SDdfno3 <- SDdfnoska %>%
  filter(county != "Stockholm")

SDdfno3def <- SDdfno3 %>%
  filter(county != "Västra Götaland")

rc31 <- polr(as.ordered(meeting_num) ~ Gender + Policyinfluence + Ambitioncand +
               Partyrecruitment + Newfriends + Stigma +
               age + agesquared + Tertiaryedu + Labourforce + yearsmember,
             data = SDdfno3def)
modelsummary(rc31, stars = TRUE, vcov = ~county, coef_rename = coefficientsnames)

rc32 <- polr(as.ordered(online_num) ~ Gender + Policyinfluence + Ambitioncand +
               Partyrecruitment + Newfriends + Stigma +
               age + agesquared + Tertiaryedu + Labourforce + yearsmember,
             data = SDdfno3def)
modelsummary(rc32, stars = TRUE, vcov = ~county, coef_rename = coefficientsnames)

rc33 <- polr(as.ordered(flyers_num) ~ Gender + Policyinfluence + Ambitioncand +
               Partyrecruitment + Newfriends + Stigma +
               age + agesquared + Tertiaryedu + Labourforce + yearsmember,
             data = SDdfno3def)
modelsummary(rc33, stars = TRUE, vcov = ~county, coef_rename = coefficientsnames)

rc34 <- glm(Partyoffice ~ Gender + Policyinfluence + Ambitioncand +
              Partyrecruitment + Newfriends + Stigma +
              age + agesquared + Tertiaryedu + Labourforce + yearsmember,
            data = SDdfno3def, family = "binomial")
modelsummary(rc34, stars = TRUE, vcov = ~county, coef_rename = coefficientsnames)
PseudoR2(rc34)

rc35 <- glm(Candidature ~ Gender + Policyinfluence + Ambitioncand +
              Partyrecruitment + Newfriends + Stigma +
              age + agesquared + Tertiaryedu + Labourforce + yearsmember,
            data = SDdfno3def, family = "binomial")
modelsummary(rc35, stars = TRUE, vcov = ~county, coef_rename = coefficientsnames)
PseudoR2(rc35)

rctable7 <- list(rc32, rc31, rc33, rc34, rc35)

modelsummary(rctable7, stars = TRUE, vcov = ~county, coef_rename = coefficientsnames,
             output = "Table RC SD - No 3.docx")


### PCA ###

rc36 <- lm(PCA_activism ~ Gender + Policyinfluence + Ambitioncand +
             Partyrecruitment + Newfriends + Stigma +
             age + agesquared + Tertiaryedu + Labourforce + yearsmember + Party,
               data = PRRdataframe2, weights = df_weights)
modelsummary(rc36, stars = TRUE, vcov = ~Party, coef_rename = coefficientsnames)

modelsummary(rc36, stars = TRUE, vcov = ~Party, coef_rename = coefficientsnames,
             output = "RC PCA.docx")
